Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2023 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
Chd|====================================================================
Chd|  RUSER36                       source/elements/spring/ruser36.F
Chd|-- called by -----------
Chd|        RFORC3                        source/elements/spring/rforc3.F
Chd|-- calls ---------------
Chd|        GET_U_FUNC                    source/user_interface/ufunc.F 
Chd|        GET_U_GEO                     source/user_interface/upidmid.F
Chd|        GET_U_MAT                     source/user_interface/upidmid.F
Chd|        GET_U_MNU                     source/user_interface/upidmid.F
Chd|        GET_U_PNU                     source/user_interface/upidmid.F
Chd|====================================================================
      SUBROUTINE RUSER36(
     1             NEL     ,IPROP   ,UVAR   ,NUVAR  ,FR_WAVE,
     2             FX      ,FY      ,FZ     ,XMOM   ,YMOM   ,
     3             ZMOM    ,E       ,OFF    ,STIFM  ,STIFR  ,
     4             VISCM   ,VISCR   ,MASS   ,XINER  ,DT     ,
     5             XL      ,VX      ,RY1    ,RZ1    ,RX     ,
     6             RY2     ,RZ2     )
C-------------------------------------------------------------------------
C  PROJET PREDIT - MECALOG  (NOVEMBRE98 - NOVEMBRE 2000)
C 
C  ROUTINE USER POUR RESSORT SPECIFIQUE PREDIT
C
C  Modele : A . Combescure
C  Implementation initiale :   Q.Zeng / C. Benoit
C  Mises a jour : Q. Zeng / F. Delcroix
C-------------------------------------------------------------------------
C   I M P L I C I T   T Y P E S
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
C----------------------------------------------------------
C   D U M M Y   A R G U M E N T S   A N D   F U N C T I O N
C----------------------------------------------------------
      INTEGER NEL,NUVAR,IPROP,
     .        GET_U_PNU,GET_U_PID,GET_U_MID,GET_U_MNU,
     .        KFUNC,KMAT,KPROP
      MY_REAL
     .   UVAR(NUVAR,*),DT ,
     .   FX(*), FY(*), FZ(*), E(*), VX(*),MASS(*) ,XINER(*),
     .   RY1(*), RZ1(*), OFF(*), XMOM(*), YMOM(*),
     .   ZMOM(*), RX(*), RY2(*), RZ2(*),XL(*),
     .   STIFM(*) ,STIFR(*) , VISCM(*) ,VISCR(*) ,FR_WAVE(*) ,
     .   GET_U_MAT, GET_U_GEO, GET_U_FUNC
      EXTERNAL GET_U_MNU,GET_U_PNU,GET_U_MID,GET_U_PID,
     .         GET_U_MAT,GET_U_GEO, GET_U_FUNC
      PARAMETER (KFUNC=29)
      PARAMETER (KMAT=31)
      PARAMETER (KPROP=33)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER I,KCST,N0,
     .        IMAT1,IPROP1,IUTYP1,
     .        IMAT2,IPROP2,IUTYP2
      INTEGER J,N,NINDX,NMAX,INDEX(MVSIZ),IFUNC
      MY_REAL 
     .        FAC,RHO,AREA,IXX,IYY,IZZ,IMYZ,YOUNG,G,
     .        AREA1,IXX1,IYY1,IZZ1,RHO1,YOUNG1,G1,
     .        AREA2,IXX2,IYY2,IZZ2,RHO2,YOUNG2,G2,
     .        K11,K22,K26,K33,K35,K44,K55,K5B,K66,K6C
      MY_REAL R_Y1,R_Y2,R_Z1,R_Z2,R_Y,R_Z,R_X,FACX,FACY,FACZ,
     .        AY1,AY2,AZ1,AZ2,AY,AZ,BY1,BY2,BZ1,BZ2,BY,BZ,G3,
     .        CX1,CX2,CX,AREA_Y(MVSIZ),AREA_Z(MVSIZ),
     .        YLD(MVSIZ),DRE,DRG,DPLA_I,DPLA_J(MVSIZ),DFE,DFG,
     .        ERR,F,DF,YLD_I,C1,FN,FM,DFN,DFM,Y1,Y2,
     .        PNX,PNX2,PNY,PNY2,PNZ,PNZ2,PMX,PMX2,PMY,PMY2,PMZ,PMZ2,
     .        NX2,NY2,NZ2,MX2,MY2,MZ2,TEMPY,TEMPZ,XL3,KTRAN,KROT,
     .        NX(MVSIZ),NY(MVSIZ),NZ(MVSIZ),YOUNGM(MVSIZ),GM(MVSIZ),
     .        MX(MVSIZ),MY(MVSIZ),MZ(MVSIZ),SVM(MVSIZ),H(MVSIZ),
     .        EPSVXX(MVSIZ),EPSVXY(MVSIZ),EPSVXZ(MVSIZ),
     .        EPSCBX(MVSIZ),EPSCBY(MVSIZ),EPSCBZ(MVSIZ),
     .        SIGNXX(MVSIZ),SIGNXY(MVSIZ),SIGNXZ(MVSIZ),
     .        MOMNXX(MVSIZ),MOMNYY(MVSIZ),MOMNZZ(MVSIZ),FXAP,
     .        GAMA(MVSIZ),GAMAC(MVSIZ),XL2(MVSIZ),DEVOL,FXAV,
     .        PC1,DC1,PR1,PS1,DC2,PR2,PS2,DC,PR,PS,UEQ,DTEMP,
     .        SIG01,SIG02, SIG0,HPLA,H1,H2,M1,M2,M,DSIGY,CC1,
     .        FAC1,FAC2
      DATA NMAX/10/
C-----------------------------------------------
      KCST = 0
C
      IPROP1 = GET_U_PNU(1,IPROP,KPROP) 
      AREA1  = GET_U_GEO(2,IPROP1)
      IXX1   = GET_U_GEO(3,IPROP1)
      IYY1   = GET_U_GEO(4,IPROP1)
      IZZ1   = GET_U_GEO(5,IPROP1)
      R_Y1= GET_U_GEO(6,IPROP1)
      R_Z1= GET_U_GEO(7,IPROP1)
      IPROP2 = GET_U_PNU(2,IPROP,KPROP) 
      AREA2  = GET_U_GEO(2,IPROP2)
      IXX2   = GET_U_GEO(3,IPROP2)
      IYY2   = GET_U_GEO(4,IPROP2)
      IZZ2   = GET_U_GEO(5,IPROP2)
      IMAT1  = GET_U_PNU(1,IPROP1,KMAT) 
      IFUNC  = GET_U_MNU(1,IMAT1,KFUNC)
      YOUNG1 = GET_U_MAT(7,IMAT1)
      G1     = GET_U_MAT(6,IMAT1)
      RHO1   = GET_U_MAT(0,IMAT1)
      AY1    = GET_U_MAT(8,IMAT1)
      AZ1    = GET_U_MAT(9,IMAT1)
      BY1    = GET_U_MAT(10,IMAT1)
      BZ1    = GET_U_MAT(11,IMAT1)
      CX1    = GET_U_MAT(12,IMAT1)
      DC1    = GET_U_MAT(13,IMAT1)
      PR1    = GET_U_MAT(14,IMAT1)
      PS1    = GET_U_MAT(15,IMAT1)
      SIG01  = GET_U_MAT(16,IMAT1)
      H1     = GET_U_MAT(17,IMAT1)
      M1     = GET_U_MAT(18,IMAT1)
      FAC1   = GET_U_MAT(19,IMAT1)
      IMAT2  = GET_U_PNU(1,IPROP2,KMAT)
      YOUNG2 = GET_U_MAT(7,IMAT2)
      G2     = GET_U_MAT(6,IMAT2)
      RHO2   = GET_U_MAT(0,IMAT2)
      AY2    = GET_U_MAT(8,IMAT2)
      AZ2    = GET_U_MAT(9,IMAT2)
      BY2    = GET_U_MAT(10,IMAT2)
      BZ2    = GET_U_MAT(11,IMAT2)
      CX2    = GET_U_MAT(12,IMAT2)
      DC2    = GET_U_MAT(13,IMAT2)
      PR2    = GET_U_MAT(14,IMAT2)
      PS2    = GET_U_MAT(15,IMAT2)
      SIG02  = GET_U_MAT(16,IMAT2)
      H2     = GET_U_MAT(17,IMAT2)
      M2     = GET_U_MAT(18,IMAT2)
      R_Y2   = GET_U_GEO(6,IPROP2)
      R_Z2   = GET_U_GEO(7,IPROP2)
      RHO    = HALF*(RHO1+RHO2)
      AREA   = HALF*(AREA1+AREA2)
      FAC    = RHO*AREA
      IXX    = HALF*(IXX1+IXX2)
      IYY    = HALF*(IYY1+IYY2)
      IZZ    = HALF*(IZZ1+IZZ2)
      R_Y = HALF*(R_Y1+R_Y2)
      R_Z = HALF*(R_Z1+R_Z2)
      R_X = HALF*(R_Y+R_Z)
      CX  = HALF*(CX1+CX2)
      AY  = HALF*(AY1+AY2)
      AZ  = HALF*(AZ1+AZ2)
      BY  = HALF*(BY1+BY2)
      BZ  = HALF*(BZ1+BZ2)     
      IMYZ   = MAX(IYY,IZZ)
      YOUNG  = HALF*(YOUNG1+YOUNG2)
      G      = HALF*(G1+G2)      
      FACX   = R_X/MAX(IXX,EM20)
      FACY   = R_Y/MAX(IYY,EM20)
      FACZ   = R_Z/MAX(IZZ,EM20)
      TEMPY  = G*AREA/MAX(TWELVE*YOUNG*IYY,EM20)
      TEMPZ  = G*AREA/MAX(TWELVE*YOUNG*IZZ,EM20)
      G3     = THREE*G
      N0     = 11
      DC     = HALF*(DC1+DC2)
      PR     = HALF*(PR1+PR2)
      PS     = HALF*(PS1+PS2)
      SIG0   = HALF*(SIG01+SIG02)
      HPLA   = HALF*(H1+H2)
      M     = HALF*(M1+M2)     
C
C----  UVAR(N0,I) -> PLA , N0+1 -> SY0 , N0+2 -> EXX_P , N0+3 -> EYY_P
C----  UVAR(N0+4,I) -> ENDOMMAGEMENT
C
      DO I=1,NEL
        YOUNGM(I)  = YOUNG*(1-UVAR(N0+4,I))
        GM(I)      = G*(1-UVAR(N0+4,I))
C== ON CALCULE D'ABORD LES VITESSES DES DEFORMATIONS GENERALISEES ====
        DTEMP = ONE/MAX(XL(I),EM20)
        EPSVXX(I) = VX(I)*DTEMP
        EPSVXY(I) = -HALF*(RZ1(I) + RZ2(I))
        EPSVXZ(I) = HALF*(RY1(I) + RY2(I))	
        EPSCBX(I) = R_X*RX(I)*DTEMP
        EPSCBY(I) = R_Y*(RY2(I) - RY1(I))*DTEMP
        EPSCBZ(I) = R_Z*(RZ2(I) - RZ1(I))*DTEMP
C== ON CALCULE CONTRAINTES GENERALISEES ====
        XL2(I)    =  XL(I)*XL(I)
        AREA_Y(I) = AREA/(ONE +TEMPY*XL2(I))    
        AREA_Z(I) = AREA/(ONE +TEMPZ*XL2(I))
        SIGNXX(I) = FX(I)/MAX(AREA,EM20)+YOUNGM(I)*EPSVXX(I)*DT
        SIGNXY(I) = FY(I)/MAX(AREA_Y(I),EM20)+GM(I)*EPSVXY(I)*DT
        SIGNXZ(I) = FZ(I)/MAX(AREA_Z(I),EM20)+GM(I)*EPSVXZ(I)*DT
        MOMNXX(I) = XMOM(I)*FACX+GM(I)*EPSCBX(I)*DT
        MOMNYY(I) = YMOM(I)*FACY+YOUNGM(I)*EPSCBY(I)*DT
        MOMNZZ(I) = ZMOM(I)*FACZ+YOUNGM(I)*EPSCBZ(I)*DT
      ENDDO
C
C--- UVAR(11,I)-> DEF PLASTIQUE EQUIVALENTE -
C--- UVAR(N0,I)-> PLA, N0+1->SY0; N0+2->EXX_P; N0+3->EYY_P;
C--- UVAR(N0+3,I)-> COURBURE PLASTIQUE Y,Z UVAR(N0+2,I)-> COURBURE PLASTIQUE X -
C
      DO I=1,NEL
        IF (IFUNC /= 0) THEN
          Y1= FAC1*GET_U_FUNC(IFUNC,UVAR(N0,I),Y2)
          Y2= FAC1*Y2   
          YLD(I) = MAX(Y1,EM20)                
          M=1
        ELSE
          Y1 = SIG0+HPLA*(UVAR(N0,I))**M
          Y2 = HPLA
        ENDIF
        YLD(I) = MAX(Y1,EM20)
        H(I)   = MAX(Y2,ZERO)
        UVAR(N0+1,I)=MIN(YLD(I),UVAR(N0+1,I))
C
        PC1 = ABS(SIX*YOUNGM(I)*(UVAR(N0+2,I)/UVAR(N0+1,I)))
        C1=ONE-FOURTH*EXP(-PC1)
        GAMAC(I)=CX*THREE_OVER_4/MAX(C1,EM20)
C
        PC1 = SIX*YOUNGM(I)*(UVAR(N0+3,I)/UVAR(N0+1,I))
        CC1 =  (THREE*PI/SIXTEEN)**2
        C1=ONE - (ONE - CC1)*EXP(-PC1)
        GAMA(I)=CC1/MAX(C1,EM20)	
C
        NX(I)=SIGNXX(I)
        NY(I)=AY*SIGNXY(I)
        NZ(I)=AZ*SIGNXZ(I)
        MX(I)=GAMAC(I)*MOMNXX(I)
        MY(I)=BY*GAMA(I)*MOMNYY(I)
        MZ(I)=BZ*GAMA(I)*MOMNZZ(I)
        SVM(I)=NX(I)*NX(I)+THREE*(NY(I)*NY(I)+NZ(I)*NZ(I)+MX(I)*MX(I))+
     1         MY(I)*MY(I)+MZ(I)*MZ(I)
        SVM(I)=SQRT(SVM(I)) 
      ENDDO
        FAC2 = G/MAX(EM20,YOUNG)
      DO  I=1,NEL
C
C       TIME STEP
C
        DTEMP    = ONE/MAX(XL(I),EM20)
        KTRAN    = MAX(AREA,FAC2*AREA_Y(I),FAC2*AREA_Z(I))*DTEMP
        KROT     = FOUR *MAX(IYY*DTEMP,IZZ*DTEMP)
        STIFM(I) = YOUNGM(I)*KTRAN
        STIFR(I) = MAX( GM(I)*IXX*DTEMP,YOUNGM(I)*KROT)
        VISCM(I) = ZERO
        VISCR(I) = ZERO
        MASS(I)   = XL(I)*FAC
        XINER(I)  = XL(I)*RHO*MAX(IXX,IMYZ+AREA*XL2(I)/TWELVE)
      ENDDO
C
      NINDX=0
      DO I=1,NEL
        IF (SVM(I) > YLD(I) .AND. OFF(I) == ONE) THEN
          NINDX=NINDX+1
          INDEX(NINDX)=I
        ENDIF
      ENDDO
C
      IF (NINDX /= 0) THEN
        DO J=1,NINDX
          I=INDEX(J)
          DPLA_J(I)=(SVM(I)-YLD(I))/(G3+H(I))
        ENDDO
C
C ATTENTION, LA VERSION DU MODELE DE RESSORT (Q.ZENG/A.COMBESCURE/C.BENOIT)
C NE CONSIDERE PAS D'ERREUR CUMULEE SUR PLUSIEURS ELEMENTS...
C
C --> EXISTE VERSION AVEC TEST MULTI-ELEMENTS (ERR_MULTI/RUSER29.F)
C
C
C
        N=0
        ERR =2*EM4
        DO WHILE (ERR > EM4 .AND. N < NMAX)
          DO J=1,NINDX
            I=INDEX(J)
            DPLA_I=DPLA_J(I)
            IF (IFUNC /= 0) THEN
              YLD_I =YLD(I)+H(I)*DPLA_I
	    ELSEIF (UVAR(N0,I) > EM20) THEN
              YLD_I =YLD(I)+H(I)*M*DPLA_I*(UVAR(N0,I)+DPLA_I)**(M-1)
	    ELSE
              YLD_I =YLD(I)
	    ENDIF
            DRE =YOUNGM(I)*DPLA_I/YLD_I
            DRG =G*DPLA_I/YLD_I
            PNX  =ONE/(ONE+DRE)
            PNY  =ONE/(ONE+THREE*DRG*AY**2)
            PNZ  =ONE/(ONE+THREE*DRG*AZ**2)
            PMX  =ONE/(ONE+THREE*DRG*GAMAC(I)**2)
            PMY  =ONE/(ONE+DRE*(BY*GAMA(I))**2)
            PMZ  =ONE/(ONE+DRE*(BZ*GAMA(I))**2)	  
            PNX2=PNX*PNX
            PNY2=PNY*PNY
            PNZ2=PNZ*PNZ
            PMX2=PMX*PMX
            PMY2=PMY*PMY
            PMZ2=PMZ*PMZ
            NX2=NX(I)*NX(I)
            NY2=NY(I)*NY(I)
            NZ2=NZ(I)*NZ(I)
            MX2=MX(I)*MX(I)
            MY2=MY(I)*MY(I)
            MZ2=MZ(I)*MZ(I)
            FN=NX2*PNX2+3.*(NY2*PNY2+NZ2*PNZ2)
            FM=3.*MX2*PMX2+MY2*PMY2+MZ2*PMZ2
            F = FN+FM-YLD_I*YLD_I
            DFE = NX2*PNX2*PNX+MY2*PMY2*PMY*(BY*GAMA(I))**2
     .                    +MZ2*PMZ2*PMZ*(BZ*GAMA(I))**2
            DFG = NY2*PNY2*PNY*AY**2+NZ2*PNZ2*PNZ*AZ**2
     .                    +MX2*PMX2*PMX*GAMAC(I)**2
            IF (IFUNC /= 0) THEN
              DSIGY = -H(I)*YLD_I*TWO
	    ELSEIF (UVAR(N0,I)>EM20) THEN 
              DSIGY = -(H(I)*YLD_I*M*(DPLA_I+UVAR(N0,I))**(M-1))*TWO
            ELSE
              DSIGY = ZERO
	    ENDIF
            DF =(-DFE*(YOUNGM(I)-DRE*H(I))/YLD_I
     .           -NINE*DFG*(G-DRG*H(I))/YLD_I)*TWO
            DF=DF+DSIGY
            ERR=ABS(F/DF)
            DPLA_J(I)=MAX(ZERO,DPLA_I-F/DF)
          ENDDO ! DO J=1,NINDX
          N=N+1
        ENDDO ! DO WHILE (ERR > EM4 .AND. N < NMAX)
C
C     CONTRAINTES PLASTIQUEMENT ADMISSIBLES
C
        DO J=1,NINDX
          I=INDEX(J)
          UVAR(N0,I) = UVAR(N0,I) + DPLA_J(I)
          DPLA_I=DPLA_J(I)
C INCREMENT N+1
          IF (IFUNC /= 0) THEN
            YLD_I =YLD(I)+H(I)*DPLA_I
          ELSEIF (UVAR(N0,I) > EM20) THEN
            YLD_I =YLD(I)+H(I)*M*DPLA_I*(UVAR(N0,I))**(M-1)
          ELSE
            YLD_I =YLD(I)
          ENDIF
          C1    = DPLA_I/YLD_I
          DRE =YOUNGM(I)*C1
          DRG =G*C1
          PNX  =ONE/(ONE+DRE)
          PNY  =ONE/(ONE+THREE*DRG*AY**2)
          PNZ  =ONE/(ONE+THREE*DRG*AZ**2)
          PMX  =ONE/(ONE+THREE*DRG*GAMAC(I)**2)
          PMY  =ONE/(ONE+DRE*(BY*GAMA(I))**2)
          PMZ  =ONE/(ONE+DRE*(BZ*GAMA(I))**2)	
C
          SIGNXX(I) = SIGNXX(I)*PNX
          SIGNXY(I) = SIGNXY(I)*PNY
          SIGNXZ(I) = SIGNXZ(I)*PNZ
          MOMNXX(I) = MOMNXX(I)*PMX
          MOMNYY(I) = MOMNYY(I)*PMY
          MOMNZZ(I) = MOMNZZ(I)*PMZ
C
          UVAR(N0+2,I) = UVAR(N0+2,I)+THREE*C1*GAMAC(I)**2*ABS(MOMNXX(I))    
          UEQ =  C1*SQRT( SIGNXX(I)**2 
     1                 +  GAMA(I)**4
     1                 *(BY**4*MOMNYY(I)**2 + BZ**4*MOMNZZ(I)**2) )
          UVAR(N0+3,I) = UVAR(N0+3,I) + UEQ      
C
C ENDOMAGEMENT
C
          IF (UVAR(N0,I) > PS) THEN
             DEVOL = DC*DPLA_I/(PR - PS)
          ELSE
	     DEVOL = ZERO
          ENDIF
C
          UVAR(N0+4,I) = UVAR(N0+4,I)+DEVOL
C
          IF (UVAR(N0,I) >= PR .OR. UVAR(N0+4,I) >= DC) THEN
            OFF(I) = ZERO
            UVAR(N0+4,I) = DC
          ENDIF  
        ENDDO ! DO J=1,NINDX
      ENDIF ! IF (NINDX /= 0)
C
C  TRANSLATE FORCES AND MOMENTS :
C
      DO I=1,NEL
        FX(I) = SIGNXX(I)*AREA*OFF(I)
        FY(I) = SIGNXY(I)*AREA_Y(I)*OFF(I)
        FZ(I) = SIGNXZ(I)*AREA_Z(I)*OFF(I)
        XMOM(I) = MOMNXX(I)*OFF(I)/FACX
        YMOM(I) = MOMNYY(I)*OFF(I)/FACY
        ZMOM(I) = MOMNZZ(I)*OFF(I)/FACZ
      ENDDO
C---
 1110 FORMAT(
     & 5X,'N. . . . . . . =',I5/,
     & 5X,'ERR. . . . . . =',E12.4/,
     & 5X,'DPPLAS . . . . =',E12.4/,
     & 5X,'F PLAS . . . . =',E12.4/,
     & 5X,'DF PLAS. . . . =',E12.4///)
 1120 FORMAT(
     & 5X,'MODULE YOUNG . . . . . ..=',E12.4/,
     & 5X,'PPLAS CUMULEE. . . . . . =',E12.4/,
     & 5X,'PPLAS TORSION . . . .  . =',E12.4/,
     & 5X,'PPLAS FLEXION. . . . . . =',E12.4/,
     & 5X,'ENDOMMAGEMENT . . . . .  =',E12.4///)
 1130 FORMAT(
     & 5X,'PPLAS CUMULEE. . . . . . =',E12.4/,
     & 5X,'DEVOL . . . .  . =',E12.4/,
     & 5X,'PS . . . . . ..=',E12.4/,
     & 5X,'PR . . . . .  =',E12.4///)
 1140 FORMAT(
     & 3X,E16.10,
     & 3X,F14.4)
C---
      RETURN
      END
